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Abstract 

General-purpose Monte Carlo event generators have become important tools 
in particle physics, allowing the simulation of exclusive hadronic final states. 
In this article we examine the Pythia 8 generator, in particular focusing 
on its parton-shower algorithms. Some relevant new additions to the code 
are introduced, that should allow for a better description of data. We also 
implement and compare with 2—7-3 real-emission QCD matrix elements, to 
check how well the shower algorithm fills the phase space away from the soft 
and collinear regions. A tuning of the generator to Tevatron data is performed 
for two PDF sets and the impact of first new LHC data is examined. 



b'ichard. corkc@thep.lu.se 
2 torbjorn@thcp.lu.sc 



1 Introduction 



The production of exclusive hadronic final states involves many aspects, ranging from cal- 
culable perturbative to incalculable nonperturbative physics. General-purpose Monte Carlo 
(MC) event generators remain vital tools in the simulation of such events. Here, typically, 
leading-order (LO) matrix elements are combined with parton showers and hadronisation 
models. Additionally, when the incoming beams are comprised of hadrons, parton distri- 
butions functions (PDFs) and multiple parton interactions (MPIs) are required to account 
for the (multi)parton content of the hadrons. Alternatively, a large number of special- 
ized generators concentrate on providing a more accurate perturbative description, e.g. by 
higher-order calculations. There need not be a contradiction between the two approaches; 
indeed a large amount of work has been done to combine them, in order to get the best 
possible overall description. 

Two of the most widely used general-purpose MC tools have been Pythia 6 [1] and 
HERWIG [2112] ■ These two programs both strive to describe the same physics, but differ in 
the models used to do so. One key difference between the two programs is in their shower 
algorithms; both are based on DGLAP evolution [1H6], supplemented with Sudakov form 
factors [7] , but while HERWIG uses angular ordering of emissions, Pythia 6 instead offers 
showers with either virtuality (Q 2 ) or transverse-momentum (p±) ordering. Even with all 
the differences between the two programs, it is remarkable the range of physics which both 
are able to describe. 

Despite the successes, both continue to evolve. They have been ported from Fortran 
to C++, but neither of these new versions is a simple rewrite. New features continue 
to be added, including changes to their shower algorithms. For example, Herwig++ [8] 
has a modified shower-ordering variable to be better suited for heavy particles. Other 
new shower developments include the release of the Sherpa event generator [U] , featuring a 
Catani-Seymour dipole shower ^0hT2] . and Vincia, based on the dipole-antenna picture [13] . 

The original Pythia 6 shower was virtuality-ordered. Combined with the MPI model 
[l4j, a large amount of Tevatron data can successfully be described. Later, a transverse- 
momentum-ordered shower was added in conjunction with an updated MPI model [T5IITB] . 
An advantage to this p± ordering is the implicit inclusion of coherence effects [T7], also 
present in an angular-ordered shower, but only crudely implemented in the Q 2 -ordered 
shower. Another feature is a dipole-style approach to recoils, as pioneered by the Ari- 
adne program [TT14T9] . The use of a p± evolution variable for both the shower and MPI 
frameworks also allows for interleaving, i.e. for a single downwards evolution in p±_ where, 
at any stage, either a shower emission or an MPI can take place. In Pythia 6, initial-state 
radiation (ISR) and MPIs are interleaved, while final-state radiation (FSR) is only consid- 
ered afterwards. The argument is that it is primarily MPI and ISR that compete for beam 
momentum. Also new is the inclusion of showers off all MPI subsystems, limited to just 
the hard scattering in the Q 2 -ordered showers. 

In Pythia 8, the Q 2 -ordered showers are no longer available, while the p^-ordered 
showers are updated j2U]. The main change is the additional interleaving of FSR along 
with ISR and MPI. It has been possible to obtain a good agreement with much data also 
in this new model, but some time ago it was noted that problems exist in the description 
of underlying-event (UE) physics when the model has been tuned to minimum-bias (MB) 
data J2JJ. This is exemplified in Fig. [TJ which shows the transverse region charged (a) 
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Figure 1: Rick Field's leading-jet underlying-event analysis implemented in Rivet. Results 
for Pythia 8.135, default tune, for transverse region charged (a) number and (b) J2p± 
density 

number and (b) J2p±. density for Pythia 8.135 in Rick Field's analysis of the UE activity 
as a function of the leading-jet p± [22]. Although the default tune is able to describe 
minimum-bias multiplicities, the activity in the underlying event clearly rises too quickly 
with jet p±. With this in mind, we study the new shower framework of interleaved FSR 
and ISR, and introduce modest changes to both that significantly improve agreement. 

Although the shower description is an approximation of QCD primarily intended to be 
accurate in the soft and collinear limits, it can also be used to fill the phase space away 
from these regions. One direct way to examine how well the shower does away from these 
limits is through a comparison to 2 — > 3 real-emission matrix elements. We will perform 
such a study in this article, for the improved shower framework we present. 

This should be put in context. In preparing for LHC phenomenology, there has been 
a large focus on techniques to combine higher-order matrix elements with parton showers. 
One goal is to describe many hard and widely separated jets, achieved by algorithms such 
as MLM or CKKW(-L) [23H25] . while another is to combine full NLO calculations with 
showers, with algorithms such as MC@NLO and POWHEG [26H28] . Attempts to combine 
both NLO calculations and higher-order Born-level matrix elements have also been made 
[T3 1I21IH3T] . Pythia 8 comes with matching built in for a few processes [221E3], and can 
be interfaced to external matching programs [34], but QCD multijets has not been high on 
the list. An eventual goal would be to have at least a complete built-in matching of the 
shower to the QCD 2—7-3 matrix elements for the first shower emission, e.g. in the spirit 

of [33]. 

There are at least two good reasons why the default showers should still be fixed up to 
give as good a description of the three-jet (and higher-jet) region as possible. One is that 
the CKKW-L-matching approach is based on obtaining Sudakov form factors from trial 
showers, so the more accurate the shower, the more precise the matching. The other is 
that one can foresee that time-consuming matching approaches will primarily be applied to 
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the hardest interaction of an event, whereas further MPIs will be handled to the accuracy 
provided by the default showers. 

The formulation of sound physics models is central to the development of a successful 
generator, but it is not enough. By the nonperturbative, and therefore incalculable, nature 
of much of the physics, and by the uncertainty that comes from shower approximations and 
unknown higher orders in the perturbative regime, generators inevitably come to contain 
a number of free parameters. Tuning of these parameters thus becomes a key part of 
obtaining an agreement with data. Assuming jet universality, LEP data can be used to 
tune FSR and hadronisation. The p± spectrum of Z° provides a handle on ISR (and, at 
low p±, on primordial fcjj- PDFs are a topic unto themselves, addressed further in later 
sections, still leaving the MPI model, the crosstalk between ISR and FSR, beam remnants, 
colour flow issues, diffractive physics and more. The underlying-event analyses here serve 
as vital ingredients [22j|35j|36] . Many of the tunes for the Q 2 -ordered showers of Pythia 6 
are based on Rick Field's original Tune A |37j. 

The tuning process has now been systematised to a large extent, and Rivet [3S] and 
Professor [2H] have become the tools of choice for many MC programs. For Pythia 6, 
the Professor team have released both a Q 2 -shower tune, Pro-Q20, and a p^-shower tune, 
Pro-pTO [39|I4T)]. For the p_i_-ordered showers, Peter Skands has also released the "Perugia" 
tunes, where a central parameter set is supplemented by 8 related variations that attempt 
to systematically explore different changes in the parameter sets |JT]. The current default 
tune for Pythia 8 is based on a tuning to LEP data by Hendrik Hoeth, using the Rivet + 
Professor framework, supplemented by hadron collider data comparisons by Peter Skands. 
Other tunes are often developed within the experimental collaborations (12JI23] and tunes 
that incorporate (preliminary) LHC data have also started to appear [H]. While the new 
energies probed by the LHC offer an exciting opportunity to test and constrain models 
and model parameters inside generators, data from lower energy runs (y/s = 900 GeV and 
2.36 TeV) currently appears slightly incompatible with earlier Tevatron results [45J. It 
remains to be seen if a "global tune" that can encompass both LHC and Tevatron data is 
possible. 

In Section 2, the shower algorithms of Pythia 8 are reviewed, and some developments 
made to allow an improved description of underlying-event and minimum-bias data. Other 
parts of the event generation framework are also outlined, especially those not previously 
documented and relevant for tuning exercises. In Section 3, a comparison of the shower is 
made to 2 — > 3 QCD real-emission matrix elements. In Section 4, some initial tunes of the 
modified framework are made to Tevatron data, while Section 5 contains a first look at the 
new LHC data. Finally, in Section 6, a summary and outlook is given. 

2 Event-generation framework 
2.1 Interleaved evolution 

The Pythia 8 showers are ordered in transverse momentum [15], both for ISR and for FSR. 
Also, multiparton interactions (MPI) are ordered in pj_ [14J. This allows a picture where 
MPI, ISR and FSR are interleaved in one common sequence of decreasing p± values [20J. 
This is most important for MPI and ISR, since they are in direct competition for momentum 
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from the beams, while FSR (mainly) redistributes momenta between already kicked-out 
partons. The interleaving implies that there is one combined evolution equation 



that probabilistically determines what the next step will be. Here the ISR sum runs over 
all incoming partons, two per already produced MPI, the FSR sum runs over all outgoing 
partons, and p± ma , x is the p±_ of the previous step. ISR is described by backwards evolu- 
tion [16], wherein branchings are constructed from the hard process, back to the shower 
initiators. Starting from a single hard interaction, eq. §Q can be used repeatedly to con- 
struct a complete parton-level event of arbitrary complexity. 

The decreasing p± scale can be viewed as an evolution towards increasing resolution; 
given that the event has a particular structure when activity above some p± scale is resolved, 
how might that picture change when the resolution cutoff is reduced by some infinitesimal 
dp±7 That is, let the "harder" features of the event set the pattern to which "softer" features 
have to adapt. It does not have a simple interpretation in absolute time; all the MPIs occur 
essentially simultaneously (in a simpleminded picture where the protons have been Lorentz 
contracted to pancakes), while ISR stretches backwards in time and FSR forwards in time. 
The closest would then be to view eq. ([1]) as an evolution towards increasing formation 
time. 

2.2 Dipoles vs. Feynman graphs 

Before the detailed presentation of the initial- and final-state showers, it is useful to put 
them in the context of the overall philosophy. As a simple illustration, consider the process 
qg — > qg. One of the contributing graphs is a t-channel exchange of a gluon, with one color- 
anticolour being annihilated in the process, another created, and one colour flowing through 
from the initial to final state, Fig. EK This can be transformed to Fig. [2b, illustrating how 
the partons move out from the interaction vertex. In this picture, the incoming particles are 
represented by fictitious outgoing antiparticles; from a colour point of view, these "holes" 
also represent the beam remnants left behind by the hard interactions. The three colour 
lines of the Feynman graph here transforms into the presence of three colour dipoles. By 
the nature of the endpoints, they can be classified as final-final (FF), final- initial (FI) or 
initial- initial (II). We begin by considering the first kind, but with many of the statements 
relevant for all three. 

In planar QCD, i.e. in the limit of infinitely many colours, Nc — > oo [17], the radiation 
pattern of soft gluons becomes a simple sum of the radiation from each dipole separately. 
Interference between the dipoles are suppressed by a factor 0(1/Nq) and so can often 
be neglected. When emitted gluons are not soft, however, recoil effects come into play; 
the endpoint partons of a radiating dipole obtain shifted momenta, affecting all dipoles 
they belong to. It is therefore not possible to consider emission from the different dipoles 
independently. Instead, one should impose a well-defined common order in which possible 
emissions are considered, such as p± [TTJUH]. Each new gluon emission also means that 
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Figure 2: (a) One possible Feynman graph and colour flow for the qg — > qg process, and 
(b) the resulting colour dipoles between the scattered partons and the beam remnants 



the radiating dipole is split into two, so the number of dipoles will increase as the shower 
evolves. 

While radiation should be viewed as emanating from the dipole-as-a-whole, it is nec- 
essary to impose a strategy for how the recoil of an emission should be shared between 
the two endpoints of a dipole. There is no unique prescription, but obviously the original 
dipole direction should not change more than necessary. One possible choice is to view the 
emission as being associated with one of the two endpoint partons, the radiator, while the 
other parton acts as a recoiler. If viewed in the rest frame of the dipole, the recoiler does not 
change direction by the emission, but it obtains a reduced absolute momentum such that 
the overall four-momentum of the dipole is preserved during the emission. The radiator 
and the emitted parton obtain opposite p± kicks by the branching, but the distribution in 
azimuthal angle is isotropic, at least to first approximation. All partons, both before and 
after the emission, are put on mass shell. 

By the radiator-recoiler strategy, the radiation of a dipole is split into that of its two 
endpoints. This can be done in a smooth manner, such that radiation close to either 
endpoint is associated with that parton radiating, whereas large-angle radiation is split 
between the two in reasonable proportions. In this way, it becomes possible to describe 
gluon radiation from a quark and gluon end by the splitting kernels 

1 + z 2 

-Pq^qg = C F - - - , (2) 

2 g ^ sg 2 1-z ' 1 ' 

respectively, with a continuous transition between the two for a qg dipole. 

The dipole picture can also be extended to apply to initial-state partons. Here the effect 
of radiation is that the incoming parton has to have a larger momentum than previously 
considered, in the spirit of backwards evolution. There are then four kinds of dipole ends, 
FF, FI, IF and II, classified by whether the radiator and recoiler are in the final or initial 
state. 

For ISR, however, the dipole picture is known to have shortcomings. Specifically, it is in 
direct conflict with results from traditional Feynman-diagram calculations. As an example, 
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Figure 3: (a) Colour flow (dipole) topology for qq — > Z°g, and (b) in Pythia the emission 
of a second gluon would give a recoil to the whole existing Z°g system rather than only to 
the dipole-connected first gluon 

consider the production of a Z°. For the emissions of a first gluon the two approaches would 
agree, but this emission would lead to a picture with two IF dipoles, stretched from the 
incoming q and q to the gluon, respectively, Fig. [3^l. Any further radiation would therefore 
occur from these two dipoles, and their daughters, but never again have any impact on 
the kinematics of the Z°. With Feynman diagrams, on the other hand, the Z° takes a 
recoil that is modified as further gluon emissions are considered, and typically resummation 
techniques |48j are used to sum up the effects of infinitely many gluon emissions on the 
p± spectrum of the Z°. This is a well-known problem e.g. with Catani-Seymour (CS) 
dipoles [10J. Therefore CS-based dipole showers do not necessarily apply the CS recoil 
strategy to the letter. Furthermore, when matching procedures are used, hard emissions 
are described by matrix elements and only softer radiation would be affected by the dipole 
picture. 

Given these shortcomings, the ISR in Pythia is not based on the CS-style dipole 
picture. Instead, an ISR emission is allowed to give a recoil to all partons "downstream" of 
it, i.e. that have been created by the hard process or by previously considered emissions at 
a larger p±_ scale, Fig. [3b. This applies both to the II and IF dipole ends. In the latter case 
one would still expect a memory of the colour flow topology to show up as an azimuthal 
anisotropy of radiation, studied further in later sections. For the FF and FI types, a more 
standard dipole picture is used, but we note and address an issue with the allowed emission 
region for the FI case. As it turns out, these moderate changes to the IF and FI dipole 
handling do help address the problem with tuning to Tevatron data. 

2.3 Transverse- momentum-ordered showers 

The current Pythia parton shower orders FSR emissions in terms of a p\ evolution variable, 
with an additional energy-sharing variable z in the branching. For QCD emissions the 
DGLAP evolution equations lead to the probability for the splitting of parton a — y be 

d^=^E^*W*)d*, (4) 
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where P a ^bc are the DGLAP splitting kernels. This inclusive quantity can be turned into 
an exclusive one by requiring that, for the first ("hardest") emission below the starting 
scale Pxmax' 110 emissions can have occurred at a scale larger than p\. The no-emission 
probability is the Sudakov form factor 

\ J P± P ± b,C Jz min(p' 2 ± ) 27r J 

Its introduction turns the unnormalised distribution into a normalised one, 

d^a = ^E ^f^Pa^ciz) dz S a ( P l max , pi) , (6) 
P± b,c Z7r 

i.e. with unit integral over the full phase space. In practice, a lower cutoff, p 2 ±0 , is introduced 
to keep the shower away from soft and collinear regions, which leads to a fraction of events 
with no emissions inside the allowed region. 

For initial-state radiation (ISR), using backwards evolution, a given parton b entering a 
hard scattering is unresolved into a parton a which preceded it. Here, the parton distribu- 
tion functions, reflecting the contents of the incoming hadron, must be taken into account. 
This leads to a Sudakov with the form 

*<*. Pi,.,., pd - ex P f- t'r dpi e r:T: ' «> **w xj rf^) . 

y hi t^min(p'i) 2tt xf b (x,p' ± ) J 

(7) 

where z — x/x' — Xb/x a , and a corresponding normalised distribution 

dn = ^- ^ ^fc(^) 3.^(3.^2^) d2 5 ^,Pi max ,Pi) • (8) 

Currently both the running renormalisation and factorisation shower scales, i.e. the scales 
at which a s and the PDFs are evaluated, are chosen to be p\. 

For the following, a relevant aspect is that, while p± ordering is used both for MPI, ISR 
and FSR in Pythia, the pj_ definition is slightly different in the three components. For an 
MPI the p ± is the expected one; the transverse momentum of the two scattered partons in 
a 2 —j- 2 process, defined in a Lorentz frame where the two incoming beams are back-to- 
back. To understand why the same choice is not made for ISR or FSR, consider a q — > qg 
branching, where the p± of the emitted gluon is defined with respect to the direction of the 
initial quark. The p± as a function of the gluon emission angle 9 increases up till 90°, but 
then decreases again, p± — > for 9 — > 180°. Thus, an ordering in such a p± would classify a 
~ 180° emission as collinear and occurring late in the evolution, although it would involve 
a more off-shell propagator than an emission at 90°. It could also erroneously associate a 
l/p\ divergence with the 9 — > 180° limit. Therefore it is natural to choose an evolution 
variable that does not turn over at 90°. 

To this end, consider a branching a — >■ be (e.g. q — > qg), where z is defined as the 
lightcone momentum along the axis of a that b takes. Then 

P±lc = z(l - z ) m l - (1 - z ) m l - zm l , (9) 
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and this equation can be used as inspiration to define evolution variables 

ISR : pi evol = (1 - z)Q 2 (10) 
with ml = —Q 2 and ml = ml = 0, 

FSR : pievoi = z(l _ Z )Q 2 (11) 
with m 2 a = Q 2 and ml = m 2 c = 0, 

which are monotonous functions of the virtuality Q 2 . In the actual branching kinematics, 
z in replaced by a Lorentz invariant definition while the Q 2 interpretation is retained. 

2.3.1 Initial-state showers 

For the initial-state shower, the radiating dipole is always chosen such that the recoiler 
is the incoming parton from the other side of the subcollision. Thus the whole enclosed 
systems share the recoil of a new emission, unlike e.g. a Catani-Seymour dipole scheme, as 
discussed previously. 

The z definition in eq. ( II Op is chosen to be z = ml r /ml r , where r is the recoiler, which 
allows for a straightforward bookkeeping of the x values to be used in PDFs. The actual 
p± of b and c then becomes 



/6, ■tt--.)<?-$-=d^-J^, (12) 

ar "_LevoLmax 



where Pj_evoi max = (1 — z )Qm&x = (1 — z ) 2m ar * s kmematically possible upper limit, 
not to be confused with the starting scale for the downwards evolution in p\ emV Usually 

P±evol ^ Pievol.max' SO t na t P±b,c ~ P±evol- 

To understand how the two dipole ends combine to produce the radiation pattern, again 
consider q(l)q(2) — > g(3)Z°(4), cf. Fig. [3^. Defining 9 as the emission angle of the gluon in 
the rest frame of the process, and with the the Mandelstam variables t,u — — s(l=pcos^)/2, 
for the gluon being emitted by the incoming q, one has 

^ = = d(C ° S§ \ = (1 + cos§) d(c ° Sg \ , (13) 

Plevol Q 2 t 1-COS^ 1-COS 2 

while an emission by the incoming q is exactly the same, except for a sign flip from 
It then follows that the quark is responsible for a fraction (1 + cos 9)/ 2 of the total amount 
of radiation at an angle 9, i.e. there is a smooth transition in the middle of the event. One 
simplifying factor is that z = M%/s is independent of the radiator side. Another is that 
we have considered the process for fixed incoming partons, so that PDF weights do not 
enter. If viewed instead as backwards evolution from qq — > Z°, a different PDF ratio would 
enter for emissions on the two sides, but the rest would be the same. The simple relations 
between the evolution and the Mandelstam variables makes this shower convenient for some 
matching tasks [32J . Further details of the kinematics, also for massive quarks, can be found 
in p2]. 

The direct relation between p_i_ ev oi an d emission angle (for fixed z) implies that typically 
central emissions will be considered first, and then successively emissions closer in angle or, 
equivalently, rapidity to the two beam directions. This ordering need not be strict however, 
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since z will also vary from one emission to the next. This implies that the resulting dipoles 
can zigzag back and forth in rapidity. We note that in CCFM-based [191150] approaches [5T] 
(as well as in angular-ordered showers, of course) emissions are ordered in rapidity, but 
then need not be ordered in p±. Hard emissions tend to be ordered in both, however, 
and it is only in the soft region that a random walk in p± sets in. While the amount of 
emitted partons need not be incorrect in our p^-ordered approach, there is a worry that 
the zigzagging colour connections will give too large dipoles, which then gives more FSR 
and hadronisation activity than if colours were angularly ordered. Ultimately this is an 
issue that could be studied by the rate of forward jets and forward particle production in 
general. Different modifications to the existing algorithm could be envisioned. For now, 
however, we will content ourselves with an option where rapidity-unordered ISR emissions 
are vetoed, subsequent to the first emission off each dipole end. That is, while p± is always 
ordered, the rapidity may or may not be. 

Perturbation theory breaks down in the p± — » limit, and thus some cutoff procedure 
is necessary for showers. The simplest is a fixed lower p± m i n cutoff. However, a more 
attractive possibility is a smooth damping. Actually, the damping introduced within the 
MPI framework, see below in Sec. \'2A.2\ could be viewed as a reduced "resolved" partonic 
content of the two incoming hadrons, and thus be the same for the ISR description. If so, 
the shower rate should be reduced roughly like 

dp 2 ± p\ P 2 ±o+pI 
We will allow for either possibility in the following. 



2.3.2 Azimuthal asymmetries 

For an ISR emission, as described previously, the dipole is always taken as stretching 
between the two incoming partons. In the colour-dipole picture, however, it will involve a 
colour line that is either stretched along the beam (collision) axis, or out to a final-state 
parton. This can lead to azimuthal asymmetries, as follows. In the first emission from a 
qq — > Z° event, the colour line stretched between the two beams means there is no preferred 
azimuthal direction, and emissions are isotropic in ip. For a second emission, the colour 
line is instead now stretched from the radiating incoming parton to the final state. This 
colour dipole configuration leads to an enhancement of radiation near the azimuthal angle 
of this final-state parton. For the emission of a soft gluon from this dipole, and with a 
convenient separation of radiation between the two dipole ends, it is possible to derive 
what this distribution should look like [52]. To be specific, if the radiator is moving along 
the +z direction and its colour partner along a (reasonably small) polar angle # par and 
azimuthal angle <£> par , then soft emissions should be distributed roughly like 

dV 1 — rcos(Aip) . . 

T~ K T^~2 — o Ta V wlth Al P = l P- Wr , (15) 

dtp 1 + r z — 2r cos(A<^J 

where r = 6/6 par . For small r the bias is tiny, but then increases as r — > 1. It turns partly 
negative for r > 1, and there integrates to — the angular-ordering condition. It must 
be stressed, however, that the behaviour near r — >■ 1 depends on how the dipole is split 
into two ends, and also obtains corrections when the energy of the radiated parton is not 
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vanishingly small. As already discussed, we also optionally allow for cases where emissions 
are not rapidity-ordered, meaning r > 1, where azimuthal anisotropies again should fall off 
in the limit of a large r. Instead, we define an alternative 

1 ^par 

where N is a free parameter to be tuned later. It has a natural value of 1/2 by the above 
small-angle behaviour, but we allow for somewhat different values, though never getting 
close to unity. 

Another reason for anisotropies is the plane polarization of the gluon. That is, the 
planes defined by the production and by the decay of a gluon, respectively, tend to be 
correlated [52] . If the gluon is produced in a g — > gg branching then a g — > gg decay leads 
to a small tendency for the two planes to line up. By contrast, a g — > qq branching gives a 
strong tendency for the two planes to be at 90° to each other. The branching g — > qq is less 
frequent than g — > gg, so to first approximation it is expected that the net effect should 
be vanishingly small [53] . This is also what we observe, both when applied for ISR and for 
FSR. 



2.3.3 Final-state showers 

For FSR, the dipole picture is again used but, unlike ISR, the recoil of emissions is here 
assumed to be taken entirely by the other end of the colour dipole, as already discussed in 
Sec. 12.21 For the case of FF dipoles all relevant kinematics equations are already published 
[15], and we only sketch some relevant points here. 

In the rest frame of a dipole ar, the branching a — > be at an evolution scale p_i_ ev oi 
implies that parton a acquires a virtuality Q 2 = p 2 ± evo \/(z(l — z)). Its energy is thereby 
increased from the original m ar /2 to {m 2 ar + Q 2 )/(2m ar ), and the recoiler energy reduced 
accordingly. Parton b then takes a fraction z of the increased radiator energy, and c takes 
1 — z. This choice gives an exact match e.g. to the singularity structure of the well-known 
qq — > q(l)q(2)g(3) branching in e + e~ annihilation: 

d Picvoi, q dz q + dpj cvoli g dzq _ dxi dx 2 + dx 1 dx 2 _ dx 1 dx 2 ^ 

P±evol,q l-*q Plcvol.q 1 ~ ^ (1 ~ X 2 )x 3 (1 - X X )x Z (1 - X X ) (1 - X 2 ) 

with %i = 2Ei/E CM . In the soft-gluon limit, x 3 — > 0, the kinematics simplifies to 1 — x 2 oc 
m^g oc 1 — cos 9 and 1 — x\ oc m| g oc 1 + cos 9, where 9 is the emission angle of the gluon with 
respect to the original qq axis. Thus one obtains exactly the same continuous transition 
from q- to q-emitted gluons as for ISR, eq. (fl3|) . For nonvanishing gluon energy, the angular 
dependence becomes more complicated, but retains the same qualitative features. 

The relation between the formal p± evo i and the p^_ b = p± c of the daughters with respect to 
the mother direction is somewhat lengthy. It shares the same physical properties as eq. (1121) 
for ISR; at small values p\ em \ and p\ bc follow suit, and so correspond to identical l/p\ 
singular behaviours, but the latter then turns around and vanishes when p^evoi approaches 
the kinematical limit. 

While FF dipoles are familiar from and constrained by e + e~ physics, the FI ones are 
less well studied. In the Pythia 6 shower, FSR was deferred until after all ISR, such that 
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the FI dipoles at the end of the ISR evolution were typically characterised by a low pj_ 
scale, meaning they did not radiate much. Indeed, it was shown that disregarding emission 
from them altogether had only a tiny influence on event properties. Instead, with ISR and 
FSR interleaved, typically each event contains one or a few FI dipole ends that should be 
evolved all the way from the p± scale of the hard interaction. 

The kinematics of an FI branching gives some differences relative to an FF one. In the 
dipole rest frame a fraction Q 2 /m 2 ar of the recoiler energy is given from the recoiler to the 
emitter, as above. That is, the emitter four-momentum is increased to 

Q 2 

Pa=Pa + -^TPt , (18) 
m lr 

which is valid in any Lorentz frame. But the recoiler is not a final-state particle, so the 
increase of a momentum is not compensated anywhere in the final state. Instead the incom- 
ing parton that the recoiler represents must have its momentum increased, not decreased, 
by the same amount as the emitter. That is, its x value needs to be scaled up as 

K = (l + J) x r . (19) 

The dipole mass m ar and the squared subcollision mass s are increased in the process, the 
latter by the same factor as x r . As with ISR, the increased x value leads to an extra PDF 
weight 

x r f r (x r , Pjj (20) 
x r f r (x r , Pj_) 

in the emission probability and Sudakov form factor. This ensures a proper damping of 
radiation in the 1 limit. 

The above equations are valid in the massless limit. Incoming particles are always 
bookkept as exactly massless, so that p' r remains parallel with the beam axis if p r was. If 
the recoiler has an on-shell mass m 2 , then its energy is already from the onset bigger than 
the recoiler one in the dipole rest frame, and the momentum and x are not rescaled by as 
much 

< = ( 1 + Q l~ m \ ) x r . (21) 

\ m ar- m a) 

The issue that remains is how the emission off the FI dipole end should be damped as 
a function of location along the dipole. It has already been shown that there is a natural 
fall-off at around 90° in the rest frame of the dipole. For an FF and II dipole this offers 
a convenient subdivision, with the two ends together providing a smooth coverage. What 
must now be asked is how to subdivide FSR and ISR in an FI dipole. The natural choice 
is to view the dipole in a boosted+rotated frame where the radiator is along the +z axis, 
with an energy equal to its original p±, and the recoiler along —z with energy m 2 ar /(Apx). 
This is close to the Breit frame definition in Deeply Inelastic Scattering, with p± as the 
hard scale instead of (the DIS variable) Q/2 and m 2 r instead of W 2 . The Breit frame offers 
a convenient separation between asymmetric current and target hemispheres. 

A more precise motivation is offered by the fractal picture of dipole emission. While 
originally formulated for FSR [51], the same kinematics argument can be applied here, 
as follows. Consider a dipole stretched between two partons i and j along the ±z axis. 
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At some arbitrary (e.g. cutoff) scale m the rapidity range open for emissions is then 
Ay = ln(m^-/mo). Assume a gluon k is emitted at an angle 6, sufficiently soft that recoil 
effects can be neglected. Then the new range(s) allowed subsequently is 

a,' = m«/ m g) + HKIO - h ( 2E ' m m ; mS > 2EsMl J "**> ) 

, (AE l E j E 2 sm 2 0)\ nl . . . 

= In — y k , )=Ay + 2\n{p Lk /m ). (22) 

Here the unchanged Ay term can be associated with radiation essentially along the z axis. 
In our case this would correspond not only to II dipoles but also to the IF ends of FI dipoles, 
which are both considered as a continuous chain along the whole allowed rapidity range. 
The last term is the additional radiation to be associated with the two new FI dipole ends. 

In summary, what we see is that the p± scale of the hard emission serves a double 
purpose. On the one hand, it sets the maximum scale for subsequent emissions, on the 
other it delimits the phase space that emissions from an FI dipole end can populate without 
doublecounting. 

This clarification especially concerns cases when the two relevant scales, the dipole mass 
m ar and the radiator scale p± are well separated. Actually, this is precisely what you expect 
to happen by t-channel gluon exchange; a (relatively) small change of momentum but a 
complete swap of colours. So in a process such as ud — > ud the outgoing u will be hooked 
up with the beam remnant of the incoming d, which most of the time will be in the other 
hemisphere of the event, with m ar ^> p±. Such a dipole would have a large rapidity range 
to radiate inside, if hemispheres are separated in the dipole rest frame, and produce a 
doublecounting of FSR and ISR effects. As we will show, indeed it does appear to make a 
difference in event properties. 

Finally, we need to specify a convenient expression for the damping at the border be- 
tween the two uneven hemispheres. Again consider the "Breit" -frame arrangement, with 
parton a having a momentum p± along the +z axis. In the branching a — > be parton b 
takes the fraction z, where we assume z < 1/2, so that b is the softer one. Both have 
transverse momentum p±b tC . If z is interpreted as lightcone momentum then rapidity for 
b is seen to correspond to p±b, c = P+b = zp +a = z2pj_. The complete kinematics is more 
messy, but qualitatively it holds that pj_b, c — z(l — z)p± is a reasonable measure of rapidity 
0, give or take factors of 2, and with 1 — z inserted for symmetry. We therefore introduce 
a suppression weight 

w = z ( l - Z )P± = P±P±b,c = P±P±b,c , 23 , 

z{\ - z)p_i_ + p ±btC p±p ±b , c + p ±b J (z(l - z)) p±p± b ,c + Q 2 

to smoothly kill the unwanted radiation in the backward hemisphere. Note that this 
damping comes on top of the "natural" subdivision of radiation in the dipole rest frame. 
Whichever one is more restrictive will dominate. 

2.4 Other model details 

In later sections, parameters relating to other parts of the event generation framework 
enter when tuning to data. In the rest of this section, then, these parts of the generator 
are outlined. 
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2.4.1 Total cross sections 



The total cross section for pp collisions (excluding the Coulomb contribution) is subdivided 
into elastic and inelastic ones, with the latter further subdivided into diffractive and non- 
diffractive contributions. Many diffractive topologies are allowed, but the most common 
and the only ones currently implemented in Pythia are single and double diffraction. We 
therefore write 



with s the squared center-of-mass energy. While the class of elastic events is unambiguously 
defined, the other three are less easily distinguished. Single and double diffraction should 
be characterised by the presence of a rapidity gap, with either one or both of the incoming 
protons excited to a higher-mass system, whereas non-diffractive events should have no 
rapidity gaps. In practice, however, generated diffractive systems can fluctuate in mass up 
to the kinematical limit, in which case there is no gap, while generated non-diffractive events 
can still display gaps, through fluctuations in the hadronisation process. The experimental 
definition of diffractive events therefore does not have a one-to-one correspondence with the 
underlying picture in a generator. This opens up the possibility for alternative experimental 
definitions, but it is important that the exact one is clearly specified. 

The "minimum bias" subset is also poorly defined as all events that are triggered within 
the context of some specific detector setup. In practice, elastic and low-mass diffractive 
events are excluded, leaving inelastic non-diffractive events with a modest high-mass diffrac- 
tive contamination. Within the Pythia generator "minbias" is used as a convenient short- 
hand for the more clumsy and nondescriptive "inelastic non-diffractive" classification. 

The different cross sections can be calculated from Regge theory, once a few free param- 
eters have been extracted from data itself. This works fine at lower, fixed-target energies. 
It also appears to work for the total cross section up to higher energies, as e.g. in the DL 
parameterisation used by us [55]. In this approach the diffractive (and elastic) cross sec- 
tions increase too fast with energy, however, and ultimately exceed the parameterised total 
cross section. This signals the need for higher-order corrections to dampen the growth. 
The parameterisations developed for Pythia [SUEZ] introduce such corrections in an em- 
pirical fashion. In a program like Phojet a more sophisticated eikonalisation approach is 
chosen instead, which leads to a slower growth of diffractive cross sections than in Pythia. 
Other arguments have also been raised why the diffractive cross section should cease to 
increase [HUES] • A recent ATLAS study of the diffractive cross section at 7 TeV (although 
not currently fully corrected for detector effects) suggests that Pythia is on the high side, 
though not overly so, while Phojet rather undershoots [60J. 

It seems likely that some tuning downwards of diffractive rates will be required. While 
it has been possible for a user to set the respective cross sections by hand, this would have 
to be done separately for each energy. We have now introduced an alternative simple mech- 
anism to dampen the growth relative to the parameterisation of [56] , without (significantly) 
affecting fixed-target phenomenology. For each of the diffractive rates it is thus possible to 
specify a maximum value that will be approached asymptotically 



We do not propose that cr J mod (s) should approach asymptotia at current energies, only that 




(24) 




(25) 
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the increase should be slowed down. Since the total and elastic cross sections are not 
affected, a reduction in diffractive cross sections means an increase of the non-diffractive 
one. 



2.4.2 Multiple interactions 

Let us now study the inelastic non-diffractive event class, and focus on just the contribution 
from MPI. The probability for an interaction is then given by 

dP M pi 1 dor / rpu-i 1 da \ 

— = — - — exp - / — T^ d P± ■ 26 

dp± a nd dp± \ J P± a nd dp_L / 

Here, interactions are supposed to be generated in a decreasing sequence of p± values, 
p±i > p±2 > p±3 > . . ., starting from some maximum scale, like the p± scale of the hard 
interaction or the upper kinematical limit. The exponential is a standard "Sudakov" factor 
that restores unitarity. It encodes the probability that there is no interaction intermediate 
in p±_ between the previous (i — 1) and the current (i) one. 

The da/dp_i_ is given by the perturbative QCD 2 — ^ 2 cross section 

da r . r . _ , n . „ , n , da 



, ^2 [ dx! [ dx 2 fi(x 1 ,p 2 ± )f j (x 2 ,p'{)-r^ , (27) 
ap ± { j j j ap ± 

representing the convolution of the hard partonic cross section, da /dpi, with the two 
incoming parton densities, /, and fj. This cross section is dominated by t-channel gluon 
exchange, and diverges roughly as dp\/p\ (further steepened by the rise of the PDFs at 
small x) . To avoid this divergence, the idea of colour screening is introduced. The concept 
of a perturbative cross section is based on the assumption of free incoming states, which is 
not the case when partons are confined in colour- singlet hadrons. One therefore expects a 
colour charge to be screened by the presence of nearby anti-charges; that is, if the typical 
charge separation is d, gluons with a transverse wavelength ~ l/p± > d are no longer able 
to resolve charges individually, leading to a reduced effective coupling. This is introduced 
by reweighting the interaction cross section such that it is regularised according to 

da a 2 (p 2 ± ) a 2 s (p 2 ±0 +pl) , , 

^ —a ► ~tt2 — , ' y^) 



dpi p\ (pi +p 2 ±y 

where p±o (related to 1/d above) is now a free parameter in the model. To be more precise, 
it is the physical cross section da /dpi that needs to be regularised, i.e. the convolution 
of da /dpi with the two parton densities, eq. ( 12T|) . One is thus at liberty to associate 
the screening factor with the incoming hadrons, half for each of them, instead of with the 
interaction. Such an association also gives a recipe to regularise the ISR divergence, as 
already noted. 

Not only p± itself, as determined e.g. from Tevatron data, comes with a large un- 
certainty, but so does the energy scaling of this parameter. The ansatz for the energy 
dependence of p±o is that it scales in a similar manner to the total cross section, i.e. driven 
by an effective power related to the Pomeron intercept (55] , which in turn could be related 
to the small-x behaviour of parton densities. This leads to a scaling 



E. 



pow 
CM 



P±o(E C u) = Pfo x l^srl , (29) 
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where -E^m * s some convenient reference energy and p*[q and E^™ are parameters to be tuned 
to data. This same scaling can also be adopted when regularising the ISR divergence, with 
the possibility to set the parameters used in the shower independently of those used in the 
MPI model. 

The so-far Poissonian nature of the framework is then changed, first by requiring that 
there is at least one interaction, such that there is a physical event, and second by including 
an impact parameter, b. For a given matter distribution, p(r), the time- integrated overlap 
of the incoming hadrons during collision is given by 



after a suitable scale transformation to compensate for the boosted nature of the incoming 
hadrons. This matter profile, p(r) may then provides new model parameters within the 
MPI framework. In the original MPI model, accompanying the Q 2 -ordered showers of 
Pythia 6, the default selection was a double Gaussian 



such that a fraction (3 of the matter is contained in a radius 02, which in turn is embedded 
in a radius a\ containing the rest of the matter. Peter Skands subsequently noted that the 
inclusion of radiation off all scattering subsystems in an event provides a dynamical source 
of fluctuations, mirroring the effect of this double Gaussian matter distribution, and was 
able to achieve good agreement with data using only a single Gaussian profile, introducing 
no free parameters [41J. Both of these options are available in Pythia 8. A further option, 
intermediate in the sense that it has one free parameter, Ef™, is to have an overlap function 
for the convolution of the two incoming matter distributions of the form 



2.4.3 Parton densities 

Both hard processes, showers and multiple interactions make use of parton densities. Since 
these are not specified from first principles, different parameterisations are on the market. 
The CTEQ5L |61j PDF set, dating from 1999, is currently the default choice for both 
Pythia 6 and Pythia 8. A large range of PDF sets has always been available by linking 
to the LHAPDF library [52]. 

Recently, a selection of newer PDF sets has been made available directly in Pythia 8 
[63J, to make the Pythia program more easy to install and run standalone, with a repro- 
ducible behaviour. One also saves some computer time by not having to go via interfaces. 
These new sets include members from the CTEQ6 and MSTW2008 families [HUES], as 
well as specially modified sets designed for LO event generators, such as those from the 
CT09 and MRST* families [BS1EZ]- These modified sets have been adopted by experimental 
collaborations [45] . 

The reason for the modified sets is the following. Precision tests of QCD today nor- 
mally involve comparisons with NLO matrix elements convoluted with NLO PDFs. NLO 
expressions are not guaranteed to be positive definite, and do not have a simple probabilis- 
tic interpretation. For PDFs, specifically, it is well-known that the gluon has a tendency 




(30) 




(31) 




(32) 
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to start out negative at small x and Q 2 , and only turn positive by the QCD evolution 
towards larger Q 2 . In the MRST/MSTW sets this negativity is explicit, while it is masked 
in the CTEQ sets by picking an ansatz form that cannot be negative, but which still gives 
a very small gluon contribution. The combination of LO matrix elements and NLO PDFs 
thus gives an abnormally small rate of interactions at small x and p±, relative to either a 
pure LO or a pure NLO combination. This region may not be important for New Physics 
searches, but it is crucial for the modeling of MPI, i.e. minimum-bias and underlying-event 
physics, and also plays a role for showers. 

So if NLO PDFs are not convenient for use in LO generators, also LO PDFs have their 
problems. Usually NLO calculations give an enhanced rate relative to LO ones, the "K 
factor" is above unity, and offer a better description to data. This introduces a tension in 
LO PDF fits, where data sensitive to a specific x range attempts to pull in more of the total 
momentum, at the expense of other x ranges. The solution of the new modified PDFs is 
to allow the momentum sum rule to be broken (typically by 10 - 15%), so that the whole 
x range can be enhanced. 

These new PDFs therefore offer the hope to allow improved descriptions of data within 
a LO framework for MPI and showers. They are not necessarily a panacea, however |63j, 
and comparing results between the modified sets and more traditional ones will therefore 
be useful. Pythia 8 also allows one PDF set to be used for the hard process and another 
for MPI and showers, so that one could use NLO for the former while still having a sensible 
behaviour in the p± — > limit with the latter. 



2.4.4 Primordial k± 

When the combined MPI/ISR/FSR p± evolution has come to an end, the beam remnants 
will consist of the remaining valence content of the incoming hadrons, as well as any com- 
panion (anti) quarks to kicked-out sea quarks. These remnants must carry the remaining 
fraction of longitudinal momentum. Pythia will pick x values for each component of the 
beam remnants according to distributions such that the valence content is "harder" and 
will carry away more momentum. In the rare case that there is no remaining quark content 
in a beam, a gluon is assigned to take all the remaining momentum. 

The event is then modified to add primordial k±. Partons are expected to have a non- 
zero k± value just from Fermi motion within the incoming hadrons. A rough estimate 
based on the size of the proton gives a value of ~ 0.3 GeV, but when comparing to data, 
for instance the p± distribution of Z° at CDF, a value of ~ 2 GeV appears to be needed. 
The current solution is to decide a k± value for each initiator parton taken from a hadron 
based on a Gaussian whose width is generated according to an interpolation 

Ql Csoft + Q Chard m , , 

(Ji + (J mi + m 

2 2 

where Q is the hardness of a sub-collision (p± for a 2 — >■ 2 QCD process) and m its invariant 
mass, cr so ft and owa is a minimal and maximal value, and Qi and mi the respective 

2 2 

scale giving a value halfway between the two extremes. Beam remnants are assigned a 
separate width a remn comparable with cx so f t . The independent random selection of primordial 
k± values gives a net imbalance within each incoming beam, which is shared between all 
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initiator and remnant partons, with a reduction factor m/ (mi + m) for initiators of low- 
mass systems. With the fcj_'s of the two initiators of a system known, all the outgoing 
partons of the system can be rotated and Lorentz boosted to the relevant frame. During 
this process, the invariant mass and rapidity of all systems is maintained by appropriately 
scaling the lightcone momenta of all initiator partons. 

2.4.5 Colour reconnection 

The final step at the parton level, before hadronisation, is colour reconnection. The idea 
of colour reconnection can be motivated by noting that MPI leads to many colour strings 
that will overlap in physical space, which makes the separate identity of these strings 
questionable. Alternatively, moving from the limit of Nc — > oo to Nq = 3, it is not 
unreasonable to allow these strings to be connected differently due to a coincidence of 
colour. Adapting either of these approaches, dynamics is likely to favour reconnections 
that reduce the total string length and thereby the potential energy. 

In the old MPI framework, good agreement to CDF data is obtained if 90% of additional 
interactions produces two gluons with "nearest neighbour" colour connections [22]. More 
recently, an annealing algorithm has been used [68|l69j. again requiring a significant amount 
of reconnection to describe data. In Pythia 8, colour reconnection is currently performed 
by giving each system a probability to reconnect with a harder system 

where R is a user-tunable parameter and p±o is the same parameter as in eq. (128]) . 

With the above probability for reconnection, it is easier to reconnect low-p± systems, 
which can be viewed as them having a larger spatial extent, such that they are more likely to 
overlap with other colour strings. When a reconnection is allowed, the partons of the lower- 
Pi. systems are attached to the existing higher-p^ colour dipoles in a way that minimizes 
the total string length. 

Currently, however, all of this is only a convenient ansatz. More than that, given the 
lack of a firm theoretical basis, the need for colour reconnection has only been established 
within the context of specific models. 

2.4.6 Diffraction 

The diffractive treatment in Pythia 8 has been extended from the simple one implemented 
in Pythia 6, to share many more features with the sophisticated description of non- 
diffractive events [70]. This is possible by using the Ingelman-Schlein [7T] picture, wherein 
single diffraction is viewed as the emission of a Pomeron pseudoparticle from one incoming 
proton, leaving that proton intact but with reduced momentum, followed by the subsequent 
collision between this Pomeron and the other proton. The Pomeron is to first approximation 
to be viewed as a glueball state with the quantum numbers of the vacuum, but by QCD 
interactions it will also have a quark content. The Pomeron-proton collision can then be 
handled as a normal hadron-hadron non- diffractive event, using the full machinery of MPI, 
ISR, FSR and other aspects as already described above. 

The simple factorisation into a Pomeron flux times a Pomeron-proton collision should 
be viewed as an effective picture, i.e. as a first approximation to physics that is much more 
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complicated than that. Indeed the HERA studies show deviations from perfect factorisation 
[72] . and diffractive W and jet production at the Tevatron is reduced by a significant factor 
relative to naive expectations [73]. Nevertheless, this factorisation does allow what hopefully 
is a realistic picture for the bulk of the diffractive cross section. 

In the generation, the first step is to use the parameterised diffractive cross sections, 
Sec. 12.4.11 Once a single diffractive topology has been decided on, the second step is to 
use the Pomeron flux /p/ p (xp, t) to pick the fraction xp that the Pomeron takes out of the 
proton momentum, and the squared momentum transfer t of the emitted Pomeron, which 
is related to the scattering angle of the proton. To first approximation 

/p/pfap,*) « — exp(Bt) , (35) 

Xp 

where B is an energy-dependent slope parameter. The squared mass of the diffractive IPp 
system X is given by M\ = xps, so it follows that dxp/xp = dM x /M x = dlnM^. Since 
the rapidity range over which particle production can occur is proportional to ha.M x , it 
follows that all kinematically possible ranges are about equally probable, and by comple- 
mentarity, the same applies to the size of the rapidity gap. This is only a first approximation, 
however; in the program four alternative /p/ p (xp,t) shapes can be used, and all of them 
are somewhat more peaked towards lower masses than the ansatz above indicates. 

Note that only the relative shape of /p/ p (xp,t) is of interest to us. If the absolute 
normalization were taken at face value it would lead to a too rapidly increasing cr S( j(s), as 
already indicated in Sec. 12.4.11 Instead, we implicitly assume that the screening corrections 
appear as an energy- dependent but xp- and t-independent reduction factor. 

The third step is to set up the partonic state of the diffractive system. Here we make a 
distinction between low-mass and high-mass diffraction, not because physics is expected to 
be discontinuous, but because we need to consider such a wide range of diffractive masses 
Mx, where perturbation theory may be applied for the higher masses but not for the lower 
ones. A pragmatic dividing line is set at 10 GeV; below it everything is soft, above it the 
fraction of soft events smoothly drops towards zero. 

The soft description is simulated as a mix of two components. In one, the Pomeron is 
assumed to kick out one of the valence quarks of the proton, so that the partonic final state 
consists of a single string stretched between the kicked-out quark and the remnant diquark. 
In the other, the Pomeron instead kicks out a gluon. This gives a hairpin topology, where 
the string is stretched from a quark in the remnant to the kicked-out gluon and then back 
to the diquark of the remnant. At the lower end of the range, extending to 1.2 GeV, i.e. 
around the A mass, the former is expected to dominate, and then gradually the latter takes 
over. Above 10 GeV, the perturbative description starts to be used. 

The Pomeron PDFs are significantly worse known than those of the proton. Quite aside 
from the considerable measurement problems, the data usually probe a convolution of the 
Pomeron flux with its PDF. Therefore the two cannot be specified independently. It is 
not even guaranteed that the two factorise. Nevertheless a few alternatives are chosen, 
both toy-model ones and four HI sets. Default in Pythia 8 is the HI 2006 Fit B LO 
distribution. Note that the HI sets use an arbitrary normalisation of the flux, such that 
the PDFs do not come out normalised to unit momentum sum, but have a sum of roughly 
1/2. 

Strictly speaking the Pomeron PDF (convoluted with its flux) contributes to the overall 
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PDF of a proton. This implies that the diffractive component ought to be subtracted from 
the normal PDFs used for non-diffractive events. In reality the diffractive contribution 
comes out to be a tiny fraction of the event rate at medium-to-large p± scales, and so this 
doublecounting correction has been neglected for now. 

In the description of the MPI interaction probability eq. (12tjj) . <r n d = &nd( s ) plays an 
important role as a normalisation factor. For the handling of a diffractive system it should 
be replaced by &tot( s ), which unfortunately is not directly measurable. As with PDFs one 
reason is that the separation from the Pomeron flux factor is not unambiguous, and another 
is that we have allowed for a damping of cr^(s) away from the result of a naive convolution. 
Instead, an order-of-magnitude estimate is used, as follows. At low energies cr^(s) is in 
the ballpark 25-30 mb, and then slowly increases. By contrast cr^(s) is slightly below 
20 mb, i.e. roughly 2/3, as suggested by quark counting rules. The same number is now 
assumed for cr^t(s), by analogy between the two-gluon IP state and the two-quark 7r state. 
The larger colour charge of gluons could imply that the Pomeron couples more strongly to 
protons than pions do. By the same token the Pomeron wave function could be smaller, 
which would act in the opposite direction, so we will assume that such effects cancel. As 
a first guess, we therefore take 20 mb to be a reasonable number to use in eq. ([2T)j) . This 
assumes, however, that the Pomeron PDFs are normalised to unit momentum sum. Since 
the HI fits are normalised only to roughly half of this, cr n( j needs to be reduced accordingly, 
in order for the ratio to balance out. The actual number used therefore is 10 mb. As a 
sanity check we note that this gives an average charged multiplicity for diffractive systems 
at a given Mx that is comparable with non-diffractive pp events at the same energy. 

Double diffraction is handled in the same spirit. Specifically, the final state will 
here consist of two Pomeron-proton collision, each with properties as already described. 
What needs some more consideration is the form of the combined two-side Pomeron flux 
/ip/pxip/p^ipi, ^iP2, t), that in particular determines the masses of the two diffractive sys- 
tems. The old machinery [55] here gives a nontrivial correlation, while some of the alterna- 
tive forms have only been formulated for single diffraction and therefore here are used by 
minimal extension, which is likely to be too simpleminded. 

2.4.7 Hadronisation 

Once the partonic configuration has been specified, including its colour flow, the normal 
Lund string fragmentation machinery can be used to turn it into a set of primary hadrons 
[74] . Many of these are unstable and subsequently decay further. 

Here we assume that fragmentation parameters tuned to e + e~ data, notably from LEP, 
can be carried over unchanged for application to pp collisions. Such a "jet universality" 
assumption is by no means obvious. In e + e~ — > 7*/Z° — > qq there is only one original 
dipole that may radiate further, typically giving rise to one string that does not bend over 
so much that it can overlap with itself. By contrast, the MPI structure of pp events ensures 
that many colour strings overlap in space and time during the hadronisation process. The 
assumption that they then can be treated as completely independent of each other (once 
colour reconnection has been taken into account) may be overly optimistic. All kinds of 
collective effects could be imagined, bordering on those associated with a quark-gluon 
plasma. This could affect e.g. the fraction of strangeness production, which is known to 
be higher in heavy- ion collisions than in e + e~ ones [75] . For now, however, such potential 
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issues are left aside. 



3 Comparison to 3-jet matrix elements 

3.1 2 — > 3 real emission matrix elements 

In this section, we compare the first shower emission against 2—^3 real-emission matrix 
elements in QCD events. Eventually, we foresee a full matching of the first emission to these 
matrix elements, as is done e.g. for electroweak gauge boson production, but here we begin 
with some simple kinematic comparisons of the default shower. Even with such a matching, 
emissions after the first will still be handled by the standard shower machinery, and it is 
questionable if it would be used in all MPI subsystems rather than just the hardest process 
in an event. It is therefore important to have the best understanding of the shower as 
possible. There is quite some freedom in how to handle recoil effects, and when the shower 
has been designed to give sensible tails away from soft and collinear regions, the manner in 
which this is done grows in importance. Indeed, as we will show, in the soft and collinear 
regions, the showers give good agreement with the ME behaviour, while further out in phase 
space, the accuracy does degrade, but without exceedingly large deviations. In the rest of 
this section, some details of the 2—7-3 ME implementation are given, before comparisons are 
made, firstly to examine the azimuthal asymmetries introduced in Sec. 12.3. 2\ and secondly 
to study the rates and kinematics of shower emissions. 

The formulae for the squared matrix elements in the massless QCD approximation, 
summed and averaged over spin, are taken from JTB] , with crossing applied where necessary. 
These expressions are compact, and so easy to input, but the tradeoff is the lack of informa- 
tion for generating the colour structure of events; for the kinematic distributions we consider 
here, this is not a serious limitation. In what follows, the outgoing partons are labeled 3, 

4 and 5, and ordered such that p± 3 > p± 4 > p± 5 . The renormalisation scale is expected to 
relate to each vertex of the process, and by default is taken to be the geometric mean of the 
squared transverse masses of the three outgoing particles, $ = ^/m^_ 3 * m^_ 4 * rn\ b . The 
factorisation scale, instead, in a picture comparable to backwards evolution in ISR, would 
be related to the smaller scales in the event, and is taken to be the geometric mean of the 
two smallest squared transverse masses of the three outgoing particles, ^ = \Jm\ 4 * m\ b . 

A phase space generator has been written that is specially adapted for these processes, 
taking into account the important role of the soft and collinear singularities. The variables 
P±3, P±5, U3, 2/4 and are used to generate the phase space. p± 3 is first picked according to 
a d 2 p± 3 /p± 3 distribution and then p± 5 according to d 2 p± 5 /p±l. <p 3 and <£> 5 are picked flat, 
before p± A is finally reconstructed. All three rapidities are picked flat in their respective 
allowed range. Phase space cuts are placed on the two p±_ variables, p^™ m and p±™ m , 
and, additionally, there is an R = yj (Arj) 2 + (A(p) 2 separation cut, R^, applied to all 
three possible pairs of partons. p± A is not explicitly constrained, but must lie in the range 
specified by p± 3 and p± 5 . 

To validate these processes, the cross sections and kinematics have been compared to 
AlpGen [77H79] (excluding processes with no final state gluons) and MadEvent [80] (for 
processes with no final state gluons). For all numbers shown here, the CTEQ6L1 PDF set 
has been used, with a first order running of a s , with a s (M\y) = 0.13. In both AlpGen and 
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Cross section (/xb) 


Tevatron 


LHC 


Cuts 


Ref 


Pythia 


Ref 


Pythia 


1 


1404 


1409 


23938 


24164 


2 


2.786 


2.798 


171.58 


172.97 


3 


0.191 


0.191 


17.87 


17.92 



Table 1: Overall cross sections for 2 — > 3 partonic processes, for three sets of cuts (see text) 
at the Tevatron (pp, x/s = 1.96 TeV) and LHC (pp, x/s = 14TeV). Results from Pythia 
are compared to summed cross sections from AlpGen and MadEvent 

MadEvent, the factorisation and renormalisation scales have been set to follow the Pythia 
defaults. In all programs, the allowed incoming and outgoing flavours have been set to 
include the charm quark and below. In Tab. [TJ overall cross sections are given for three 
sets of cuts at Tevatron (pp, x/s = 1.96 TeV) and LHC (pp, x/s = 14 TeV) energies: 

1) Pl ™ = 5.0 GeV, p ± f n = 5.0 GeV, = 0.1 

2) P± f n = 50.0 GeV, p ± f n = 5.0 GeV, i?™ n = 0.1 

3) p± fi* = 50.0 GeV, p ± f n = 25.0 GeV, = 1.0 

The difference between the reference and Pythia cross sections is less than 1% for all 
three sets of cuts. In the PS, the amount of ISR is roughly 70%, 60% and 90% respectively, 
and so dominates, even when using a rather small -R™". In Tab. El the different classes of 
2 — > 3 processes are given, with indicative cross sections for each process for the first set of 
cuts at Tevatron and LHC energies. Overall, the cross section is dominated by gg and qg 
scattering, and the effect grows at LHC energies, in particular driven by the large increase 
in small- x gluons. 

3.2 Parton shower comparison 

The comparison we wish to make is between the 2—7-3 matrix elements and 2—7-2 
events with one additional PS emission (hereafter referred to as 2 — > 2 <g> PS). Overall, 
the generation must be such that the comparison is as fair as possible, but still allowing a 
reasonable efficiency, so the following settings have been applied: 

1) The CTEQ6L1 PDF set is used everywhere. 

2) The coupling strength for the hard process, ISR and FSR is set constant, a s = 0.130. 

3) QED radiation is disabled. 

4) The number of incoming and outgoing flavours possible is set equal for both the ME 
and PS (up to and including the charm quark). 

5) All other stages of the event generation framework are switched off. 



21 





Cross section (/xb) 


Class 


Processes 


Tevatron 


LHC 


HardQCD 


gg2ggg 


gg — > ggg 

OO OOO 


717.46 


15236.92 


HardQCD 


qqbar2ggg 


qq — > ggg 

11 OOO 


0.51 


2.47 


HardQCD 


qg2qgg 


qg — > qgg 


525.65 


6776.26 


HardQCD 


qq2qqgDiff 


qq' — > qq'g, qq' — > qq'g 

11 11 O 7 11 11 O 


45.47 


404.31 


HardQCD 


qq2qqgSame 


qq ->■ qqg, qq -> qqg 


8.92 


94.51 


HardQCD 


qqbar2qqbargDif f 


qq ->• q'q'g 


0.55 


2.68 


HardQCD 


qqb ar 2 qqb ar gS ame 


qq -> qqg 


16.08 


81.93 


HardQCD 


gg2qqbarg 


gg ->■ qqg 


71.22 


1302.09 


HardQCD 


qg2qqqbarDif f 


qg — >■ qq'q', qg — > qq'q' 


18.05 


197.95 


HardQCD 


qg2qqqbarSame 


qg -»■ qqq, qg -> qqq 


6.02 


66.01 



Table 2: The different classes of event implemented and their corresponding processes. 
Primes refer to different quark/ ant i- quark flavours and processes include all possible per- 
mutations of the initial or final state. Cross sections from Pythia are given for the first 
set of cuts as detailed in the text 

6) To improve efficiency (without affecting the outcome), a cut j3j_ m i n = 0.4 * p_i_™ m is 
used in generating the underlying 2 — > 2 events. 

Even with these settings, a couple of issues remain. One is the different scale at which PDF 
factors are evaluated in ISR, eq. ([8]). Specifically, both the numerator and denominator of 
this factor are evaluated at a scale p\, giving an implicit running a s when compared to the 
ME result. A quick check using typical scales and x values for the above cuts shows that 
this effect is not large, and we choose not to correct for this. 

More important is the presence of Sudakov form factors in the PS results. To generate 
inclusive 2 — > 2 ® PS distributions, which can be compared against the ME results, a 2 — > 2 
hard process is selected, then a parton shower emission allowed to create a 2 — > 2 <g> PS 
event. This event is then analysed to determine if it meets the required cuts, and if so, it 
is accepted and added to the statistics. The emission is then vetoed; Pythia will revert 
the event to its original 2—7-2 state, and then continue the shower evolution from the scale 
of the now- vetoed branching. Thus, a single 2—7-2 event, in the course of its evolution, 
can give rise to several 2—7-3 configurations. In fact, for an evolution variable t, ranging 
between and 1, the distribution becomes a Poissonian with an average Jq p(t)dt. By 
the veto algorithm, the Sudakov form factor is exactly compensated by the Poissonian 
distribution 

p(t) exp (- jf* p(f ) dt') • £ i QT p(0 dt'Y = p(t) , (36) 
where the sum runs over the n 2 — > 3 events that may have been generated at t' < t. 
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3.3 The a angle 



We begin by studying the azimuthal asymmetries introduced for ISR in Sec. 12.3.21 The 
strength of the asymmetry, tuned through the N parameter of eq. (Tl6i) , enters in all following 
studies, and one goal here is to fix its value. The question to be asked is, how strongly 
does the plane of an initial-state emission correlate with its colour partner in the final state, 
when present. The distribution can be calculated in the soft limit, but obtains corrections 
away from this region; hard emissions, especially in the case when strict rapidity ordering 
is not enforced (as for the first emission and optionally for subsequent ones) are expected 
to be constrained to a much lower degree. In both the ME calculation and in data, there is 
no clean separation between ISR and FSR, so we instead use the a angle observable, used 
by the CDF collaboration in order to study colour coherence effects in hadron colliders pJJ . 




At the hadron level, the idea is to search for a leading jet which is hard enough such 
that "soft" radiation is still hard enough to form secondary jets. The study is based on 
QCD events, where at least three jets must be present (after jet clustering with the CDF 
fixed-cone algorithm, -R cone = 0.7). Both the hardest and second hardest jets must lie in the 
region \q\ < 0.7 and the leading jet must have Ex_\ > HO GeV. The two leading jets must 
be approximately back-to-back in the transverse plane, \\<pi — — vr| < 7r/9, while the 
third jet must have E±$ > 10 GeV. These cuts should isolate a region in which the shower 
approximation is still valid, and azimuthal asymmetries should appear as a correlation 
between the directions of the second and third jets. 

Defining Ar] = r/ 3 — r] 2 and A(p = ip 3 — <p2, an R 2i = V Arj 2 + Aip 2 cut is introduced, 
1.1 < -R23 < The lower limit is placed due to clustering cone-size effects, while the upper 
limit avoids the region where \Aip\ approaches its maximal value of 7r. The a angle is then 
defined by 



A full discussion of the details of the different a-angle regions is given in (8TJ. Here, we 
restrict ourselves to azimuthal correlations only. The primary effect should be a shift of 
events towards the endpoints of the distribution, a — » ±7r/2, such that ISR is biased to sit 
in the plane of its colour-connected partner, Atp — > 0. It should be noted that a < will 
be slightly favoured due to the third jet having access to central rapidity regions, where 
phase space is larger. 

A similar study has also been done by the DO collaboration, but instead using (3 = 
tan -1 (sgn(i]2) A(f / Arj) as their observable JS2]. The overall conclusions of both of these 
studies is that parton shower MCs that implement colour coherence, either through the 
choice of evolution variable or through an angular veto are able to describe data. 

First, a comparison is made against the 2—7-3 ME, using the second set of cuts from 
the previous section. In this kinematic region, the softest jet is expected to come from 
radiation and should be relatively soft compared to the underlying 2—7-2 process, where 
the parton shower is expected to be most accurate. At the parton level, many of the cuts 
used in the full experimental analysis are no longer relevant, and so the only restriction 
placed is that 0.5 < -R23 < 1-0. The lower cut, although no longer necessary for clustering 
reasons, is included to enhance the effect of ISR by removing some contributions from FSR, 
which will be highest in the region of low i?23- The upper cut helps to isolate the region 
where the azimuthal weighting should have the biggest effect (r ~ 1, Sec. I2.3.2p . 





(37) 
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Figure 4: a-angle comparison for different strengths of asymmetry at (a) Tevatron (pp, 
y/s = 1.96 TeV) and (b) LHC (pp, y/s = 14TeV) energies 



Results are shown in Fig. |4] for different strengths of asymmetry, N, at (a) Tevatron 
(pp, y/s = 1.96 TeV) and (b) LHC (pp, y/s = 14 TeV) energies. For the PS, it is clear that 
when all initial-state azimuthal weighting has been switched off, there is an excess of events 
in the region of small |a|. When the weighting is switched on, the results are as expected; 
events are shifted towards the endpoints of the distribution, with the effect growing as N 
is increased. One unexpected feature is the agreement of the PS and ME at a = —tt/2 
when no azimuthal weighting is switched on. This suggests that ISR emitted toward the 
central region already has some natural bias towards small Atp. The extra weighting, then, 
overshoots in this region, but does have the desired effect at a = 7r/2, where the slope of 
the tail is brought up. In moving to higher energies, it should be noted that the same cuts 
are used, meaning phase space effects will affect the exact shape of the distribution and the 
balance between the positive and negative a regions. Considering only the first emission, 
a value of N = 0.45 gives a reasonable agreement across the whole a range and for both 
energies. 

This is not the whole story however; before event generation is complete, there is still 
the rest of the shower evolution and additional non-perturbative effects to be added. After- 
wards, some memory of the parton-level asymmetries should remain, which, in the a-angle 
observable, should show up in the same way as in the ME comparison. Additionally, at the 
hadron level we also study the the pseudorapidity distribution of the third jet (no explicit 



24 



(a) 




R23 separation cut is applied to this observable). It has been shown that colour coherence 
effects lead to a broadening in this distribution, as jets are more likely to be emitted at 
higher rapidity; a good agreement is also a sign that coherence effects are well modeled. 

The default tune of Pythia 8.142 is used, and only those settings related to the 
azimuthal asymmetries in ISR changed. The output of the generator is fed through the 
CDF a-angle Rivet analysis and bin-by-bin detector corrections added to the final results. 
The results are given in Fig. for (a) 1] 3 and (b) a. The width of the 7] 3 distribution is 
already well described, even without the extra azimuthal weighting. Increasing the strength 
of the asymmetry does lead to a very slight broadening of the distribution, but the largest 
difference is in the dip at 773 = 0. In general, the high-p^ jets will prefer to sit at central 
rapidities. Biasing Atp to smaller values, in combination with the implicit R separation cut 
from the jet algorithm, will shift the third jet slightly further out in rapidity. 

The results for the a-angle are similar to the ME comparison. At the smallest values, 
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there is a small overshoot, while at higher values, the agreement is improved. At a = tt/2, 
the tail does not come up to the levels of the data, but is no longer falling. Indeed, with the 
weighting that has been introduced, bringing the tail up to match the data is a hard task, 
indicating that additional correlations in rapidity may be involved in later stages of the 
event-generation process. Based on the CDF data, then, a stronger asymmetry is favoured, 
and a value N = 0.70 is adopted as the default value for all the following sections. 

3.4 Kinematic distributions 

The main freedom in adjusting the PS to better fill the phase space is given by the starting 
scale of the initial- and final-state showers. For those hard processes which do not con- 
tain any final-state particles which may be shower produced, the shower can begin at the 
kinematical limit, so as to populate the full phase space without risk of double counting. 
For the QCD processes considered here, where the hard process already contains particles 
which may be produced in the shower, double counting becomes an issue and the shower 
is instead started at the factorisation scale of the hard process. In the case of a 2 — > 2 
process, this defaults to the smaller of the squared transverse masses of the two outgoing 
particles. Neglecting parton masses, this implies that the p± of the shower jets are below 
the 2 — > 2 p± scale. This appears to be a reasonable starting point, but note that the recoil 
of an emitted parton can boost one of the original two jets to have a smaller p± than the 
emitted jet. Effects can also go in the opposite direction, where a parton emitted at the 
hard p± scale boosts both original partons to a higher pj_, opening an unfilled region. 

Results are shown using the default showers for the three sets of cuts defined previously, 
for the transverse momenta, rapidities and R separations of the three jets. The first set 
covers as large an area of phase space as possible, while the second isolates a region where 
the parton shower is expected to perform well. The final set should isolate a region where 
the three jets are all relatively hard, and well separated, exactly the region where the parton 
shower is not expected to perform well. In the plots that follow, we stress that the main 
aim is for a good qualitative agreement with the matrix elements. 

For the second set of cuts, results are shown in Fig. [6] for Tevatron energies (pp, \/s = 
1.96TeV) and in Fig. [7]for LHC energies (pp, a/s = 14TeV). While at the lower energy, the 
PS rate is slightly below that of the ME over most of the phase space, and at the higher, the 
ME is slightly above, all distributions are well reproduced. The small excesses in events at 
low p_i_4, p±5 > 25 GeV and -R45 > tt/2 are consistent with a slightly too high contribution 
from hard- and widely-separated jets, which we examine further below. 

In Fig. [5J the results for the third set of cuts at Tevatron energies are shown. In this 
region of hard, widely separated jets, the PS is not expected to perform so well. Indeed, 
here the PS rate is somewhat larger than the ME one. As previously, there is a noticeable 
excess in events at low p± 4 , suggesting configurations where a softer 2 — )• 2 process is 
shifted to meet the jet requirements through an ISR emission and its recoil. The rapidity 
distribution y 4 is also too broad. A simple check shows that there is no large contribution 
from an ISR emission becoming harder than the original two jets. That the shapes of p± 5 
and ?/5 are relatively well described suggests that there are effects arising from the rotation 
and boosts employed in the recoil handling, also giving an excess in the tail of Ai? 45 . We 
do not study this further at this time. Although the distributions are by no means perfect, 
we content ourselves with the fact that these events are in a disfavoured region of phase 
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Figure 6: Kinematic distributions for cut set (2) at Tevatron energies (pp, \fs = 1.96TeV): 
pif = 50.0 GeV, p±f n = 5.0 GeV, i?™ n = 0.1 



space and that the description gives the broad features of the shapes. 

In Fig. [9j the results for the first set of cuts are shown at Tevatron energies. This cut 
set is designed to be as inclusive as possible, representative of the MPI that are important 
for MB and UE physics. The distributions show a mixture of the features of the second 
and third cut sets as shown previously. In general, both the rates and shapes are well 
reproduced, a sign of the dominance of the soft and collinear regions of phase space. 

Some freedom to probe the ambiguities of the shower starting scale is offered by the 
possibility of introducing a constant factor, k, such that p_i_ m ax = k * Qf ac - This constant 
may be set separately for the initial- and final-state showers, and in all results shown so far 
has been set to k — 1. In Fig. (TUl the p± of the softest jet is shown for each set of cuts at 
Tevatron energies when k = 0.8 and k = 1.2, both for the initial- and final-state shower. 
Changes in k give variations both in rate and shape, and the agreement with the ME result 
is dependent on both the exact cuts used and the energy. Although scale variations away 
from k = 1 may give better agreement in certain regions of phase space and at certain 
energies, there is no indication that a systematic shift, higher or lower, would improve the 
overall PS description. 

Finally, as an example of the impact of these variations on final observables, the jet- 
jet azimuthal-angle distribution, as studied by DO, is shown in Fig. [TTJ Other than the 
different factors of k, all other settings remain at their default value. As with many other 
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Figure 7: Kinematic distributions for cut set (2) at LHC energies (pp, y/s = 14 TeV): 
p± f* = 50.0 GeV, p ± f n = 5.0 GeV, R™£ = 0.1 



observables, many different parameters in the overall framework will play a role, and these 
must then be adjusted to give the best description of the data. The shower k factors could 
be varied in the context of a more complete tuning to data, at the cost of an increase in 
the size of the parameter space, but based on the evidence above we choose to keep k — 1 
in all that follows. 



4 Tunes to Tevatron data 

In this section, having "validated" the now-modified PS framework in the previous section, 
we study the tuning prospects of Pythia 8 to Tevatron data. All changes to the framework 
outlined previously apply solely to the case of hadron collisions, so the existing LEP tune of 
the final-state shower and hadronisation can remain unchanged. However, even the subset 
of parameters relating only to hadron collisions is large, so the strategy is to pick a minimal 
set of key parameters, and try to find a region in this parameter space that describes 
key pieces of Tevatron data. One key goal is to show that a combined MB/UE tuning of 
Tevatron data is within reach for Pythia 8, as is possible in previous versions. It should 
be stressed that this is not designed to be a "complete" tuning of the generator. The search 
of this limited parameter space is guided by the principles we outline below, but there is 
still no guarantee that other regions of parameter space do not give a similar or better 
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Figure 8: Kinematic distributions for cut set (3) at Tevatron energies (pp, y/s = 1.96TeV): 
p L f* = 50.0 GeV, p±f n = 25.0 GeV, R™" = 1.0 



descriptions of the data. We hope that what we present here offers a good starting point 
for a more global Professor-type tuning, encompassing a wider set of data and parameters, 
as well as offering a general guide to those getting involved with such efforts. 

We also take this opportunity to move away from the CTEQ5L PDF set. As explained 
in Sec. I2.4.3[ a selection of newer sets is now available directly in Pythia 8, and we pick 
one standard LO set (CTEQ6L1) and one modified set (MRST LO**) with which to begin 
the tuning attempts. 

The following Rivet analyses have been used in the tuning process: 

CDF_2000_S4155203: Z pT measurement in Run IZ-> e+e~ events [83] 

p± and cross section measurement of e + e~ pairs in the region of 66 < m ee < 116 GeV. 

CDF_2001_S4751469: Field & Stuart Run I underlying event analysis [831 

The direction of the leading charged particle jet in an event is used to define three 
regions in r] — ip space. With Aip = <p — y?Leadingjet, the toward region is defined by 
\A(p\ < 60° (centered on the leading jet), the away region by \Aip\ > 120° and the 
transverse region by 60° < \A<p\ < 120°. All regions are constrained such that \rj\ < 1. 
Two datasets are used; one using a minimum bias trigger, the other using a "JET20" 
trigger, where one calorimeter tower cluster must have Et > 20 GeV. Results are 
given as a function of the transverse momentum of the leading jet, p 1 ^. 
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Figure 9: Kinematic distributions for cut set (1) at Tevatron energies (pp, \fs = 1.96TeV): 
pif = 5.0 GeV, p±f n = 5.0 GeV, i?™ n = 0.1 
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Figure 10: p± 5 for each set of cuts at Tevatron energies (pp, y/s = 1.96 TeV) with k = 0.8 
and k = 1.2 

CDF_2002_S4796047: Run I charged multiplicity measurement [55] 

Charged multiplicity measurements at yfs = 630 & 1800 GeV. The (p±)(V r c h) mea- 
surements have largely been superseded by CDF_2009_S8233977, below. 

D0_2004_S5992206: Run II jet azimuthal decorrelation analysis [HE] 

Distributions of Atp = (fi — (fi2, where <pi is the azimuthal angle of the hardest jet 
and if2 that of the second hardest jet. Results are given in 4 different bins of p± l . 

CDF_2008_LEADINGJETS: Run II underlying event in leading jet events [57] 
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Figure 11: Jet-jet azimuthal-angle distribution showing the changes when the shower start- 
ing scale is adjusted by a factor k = 0.8 and k = 1.2 



Regions are defined as in CDF_2001_S4751469, but with the two transverse regions 
additionally separated into "transMIN" and "transMAX" regions, based on which 
contains the lowest/highest number /Y^P± density, event-by-event. A "transDIF" 
region is defined by the difference between the transMAX and transMIN regions. 

CDF_2008_NOTE_9351: Run II underlying event in Drell-Yan [88] 

As CDF_2008_LEADINGJETS, but using a leading Z° reconstructed from a Drell-Yan 
lepton pair. 

CDF_2009_S8233977: Run II min bias cross-section analysis [89] 

Measurements of track p±, {p±)(N ch ) and f° r — 1-96 TeV. 

Our limited parameter space consists of the following: 

SigmaProcess : alphaSvalue 

a s (M z ) for the hard process. Directly affects the single-particle p± spectrum in the 
\ow-p± region and jet cross sections at high p±. For MRST LO**, a decrease of ~ 15% 
is expected due to similarly-sized increase in the PDF content (the QCD 2—7-2 cross 
section is quadratic both in PDFs and in a s ). 

SpaceShower : alphaSvalue 

a s (Mz) for the initial-state shower. Together with the regularisation parameters (see 
below), gives control over the amount of ISR activity. Although convenient to be set 
equal to a s of the hard process, above, this is not a requirement. 

SpaceShower :pT0Ref, SpaceShower :ecmRef, SpaceShower : ecmPow, 
SpaceShower : samePTasMI, SpaceShower :pTmin 

Parameters used to regularise the p± — > divergence in the initial-state shower 
(Sec. I2.3.ip . While the smooth dampening of eq. (fl4"|) is physically preferable to 
a step-function regularisation, as provided by SpaceShower :pTmin, there is currently 
no real constraint on these options. For now, we choose to fix the p_i_o ef to an energy 
independent value, rather than use the same running as used in the MI framework. 
Together with a s (M z ) above, these values are tuned mainly based on the p_i_(Z°) spec- 
trum (although this also contains large contributions from primordial k±) and the DO 
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jet azimuthal decorrelation analysis. While a change of p±o will give more activity 
only at small p± , s, a change to a s will affect the whole spectrum, although still not 
evenly due to its running. 

SpaceShower : rapidityOrder 

If switched on, rapidity-unordered ISR emissions, subsequent to the first emission off 
each dipole, are vetoed (Sec. I2.3.T]) . Studies suggest that having this option activated 
helps to dampen the rise of the underlying event, leading to a better agreement with 
data. 

Mult iplelnteract ions : alphaSvalue 

®s{Mz) for the MPI framework. This can be used to help control the amount of 
MPI activity, in conjunction with the regularisation parameters (see below). In the 
context of the MPI model there is no strong reason to vary this independently of 
SigmaProcess : alphaSvalue, and we choose to restrict them to be equal or very 
close to each other. As with ISR, a change to a s will affect the whole spectrum of 
MPI as opposed to just the soft region. This balance can have effects on e.g. (p±) (N ch ) 
distributions, which in turn depends heavily on colour reconnection effects. 

Mult iplelnteract ions :pTORef , Multiplelnteractions : ecmRef , 
Mult iplelnteract ions : ecmPow 

Parameters for the regularisation of the MPI framework (Sec. I2.4.2p . p±o is a key 
ingredient in getting both the charged multiplicity and the underlying-event distribu- 
tions correct. As we only consider Tevatron data here, the energy dependence is only 
constrained by the CDF charged multiplicity at y/s = 630 GeV. The reference energy 
is left unchanged at 1800 GeV. 

Multiplelnteractions :bProf ile, Multiplelnteractions : expPow 

Parameters used to select the matter profile used in the MPI framework (Sec. I2.4.2p . 
As noted previously, the introduction of showers off all scattering subsystems re- 
duces the need for a double Gaussian profile; instead we choose the overlap function, 
eq. ( 13"2"j) . giving one degree of freedom, E%g£. This parameter is important in the same 
distributions as p±o above; for charged multiplicities it affects high-multiplicity tail, 
while it also controls how quickly the underlying-event activity grows as a function 
of the transverse momentum of the leading jet in MB events. A value E?™ = 1 gives 
results somewhat similar to the default double- Gaussian matter profile, giving more 
fluctuations out in the high-multiplicity tail and a faster rise of the underlying event, 
while -E^pT = 2 reduces back to the single Gaussian case. We examine values between 
these two limits to try to find the best balance. 

BeamRemnants : reconnectRange 

The R parameter of eq. ( 1341) . This is primarily tuned to the latest measurement of 
(j)±)(N c i 1 ) from Run II Tevatron data, but should also give reasonable results for the 
lower energy runs. It should be noted that this parameter is non-linear; at R ~ 5, 
the amount of reconnection has almost completely saturated. 

Not included in the above list are those parameters related to primordial k± (Sec. 12.474"]) . 
For the hard component, early studies indicated that the p±(Z°) peak was well described 
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Parameter 


Tune 2C 


Tune 2M 


Tune 4C 


SigmaProcess : alphaSvalue 


0.135 


0.1265 


0.135 


SpaceShower : rapidityOrder 


on 


on 


on 


SpaceShower : alphaSvaluG 


0.137 


0.130 


0.137 


SpaceShower : pTORef 


2.0 


2.0 


2.0 


Multiplelnteractions : alphaSvalue 


0.135 


0.127 


0.135 


Multiplelnteractions :pT0Ref 


2.320 


2.455 


2.085 


Multiplelnteractions : ecmPow 


0.21 


0.26 


0.19 


Multiplelnteractions :bProf ile 


3 


3 


3 


Multiplelnteractions : expPow 


1.60 


1.15 


2.00 


BeamRemnants : reconnectRange 


3.0 


3.0 


1.5 


SigmaDif f ractive : dampen 


off 


off 


on 


SigmaDif f ractive :maxXB 


N/A 


N/A 


65 


SigmaDif f ractive :maxAX 


N/A 


N/A 


65 


SigmaDif f ractive :maxXX 


N/A 


N/A 


65 



Table 3: Parameters for Tunes 2C, 2M and 4C (described in Sec. [5]). Any parameter not 
shown is set to the default value of Pythia 8.142 

with its default setting (ohard = 2GeV), while for the soft component, variations showed 
little sensitivity in the mean p^ per parton. In Tab. [3] the parameters for two tunes, 2C 
(CTEQ6L1) and 2M (MRST LO**) are given. Any parameter not shown is set to the 
default value of Pythia 8.142. 

Selected plots are shown in Figs. [12] & [13] for tunes 2C and 2M respectively. The tunes 
are compared against data and Pythia 6 tunes Pro-Q20 and Perugia for (left-to-right, 
top-to-bottom): 

1) CDF_2002_S4796047: Charged multiplicity at ^/s = 630 GeV 

2) CDF_2002_S4796047: Charged multiplicity at = 1800 GeV 

3) CDF_2008_LEADINGJETS: Transverse charged particle density 

4) CDF_2008_LEADINGJETS: Transverse charged J2p± density 

5) CDF_2009_S8233977: (p±){N ch ) 

6) D0_2004_S5992206: Jet-jet azimuthal angle for 100 < p^ ax < 130 GeV 

It is notable that less colour reconnection is needed to match the {p±)(N c ^) data than 
in the old default. One contribution to this is the increased a s in the MPI framework, 
giving more activity at all p± values rather than just extra low-p^ activity from a change 
in Rlo- The good agreement in the DO jet-jet azimuthal distribution gives a good sign that 
the balance between contributions from ISR and MPI is well described, supplemented by 
p_i_(Z°) measurements, not shown. 

There is a large difference in the Ef™ parameter between the two tunes. In fact, in Tune 
2M, the multiplicities at both 630 & 1800 GeV do overshoot in the tails; the low value of 
Eexp is required to correctly describe the underlying event in the region p_i_i ead < 200 GeV. 
One key reason is that the the MRST PDF set combined with the retuned a s differs in 
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Figure 12: Selected Rivet plots for Tune 2C, 
and Perugia 
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Figure 13: Selected Rivet plots for Tune 2M, 
and Perugia 
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both shape and normalisation relative to CTEQ6L1, and when probing the amount of 
activity as a function of the leading jet, these differences become more apparent. A value 
Ee°p ~ 1-4 would give a better description of the overall multiplicity distributions, but 
with this setting, the average activity in this pj_i ea( j region is lower than for CTEQ6L1. In 
Tune 2M, we have compensated for this by allowing more impact-parameter fluctuations 
out to the high multiplicity tails, but expect regions in parameter space to exist which 
better address this balance. 

For both tunes, the transverse activity in the underlying event study no longer shows the 
large rise noted in Fig. [TJ Both tunes also show success in a combined MB/UE description 
at the Tevatron. Over all datasets, Tunes 2C and 2M never do significantly worse than 
Pro-Q20 and Perugia 0, and it is hoped that a full tuning, considering a larger parameter 
set, could improve agreement further. 

5 First LHC data 

5.1 Diffractive cross sections 

For the LHC data, we begin by studying the diffractive cross section and the simple scheme 
to dampen its growth outlined in Sec. 12.4.11 On the experimental side, we turn to a recent 
ATLAS study designed to enhance the diffractive content of minimum bias events [60J. 
The study relies on the Minimum Bias Trigger Scintillators (MBTS), covering two rings in 
pseudorapidity, 2.09 < \rj\ < 2.82 and 2.82 < \rj\ < 3.84. By selecting events with a hit on 
only one side of the MBTS, diffractive events are preferentially chosen. While only total 
cross sections are examined here, the ATLAS study also shows that track distributions are 
better described with the new high-mass diffractive framework of Sec. 12.4.61 It is noted that 
this study is not corrected for detector /reconstruction effects, but that such corrections are 
not expected to change the conclusions of the study. 

The particular quantity of interest is R ss = N BB /N any , where N ss is the number of events 
with a hit on exactly one side of the MBTS and N any the number with a hit on either 
side. We follow the approach of the ATLAS study, in using fixed acceptance values for the 
different event classes while varying the contributions of the different diffractive modes. 

From data, a value R ss = [4.52 ±0.02(stat) ±0.61(syst)]% is preferred, while the default 
settings from Pythia 8 gives a value R ss = 5.11%, lying just inside the error range. 
Assuming that the single- and double-diffractive components saturate at the same value, 
a damping cr^ x = <J^^ = 65 mb brings i? ss into closer agreement with the data. With 
this value, in this simpleminded model, saturation would then still be far away. Tab. H] 
gives the non-, single- and double-diffractive cross sections for different energies in both 
the default and damped scenarios. R ss values are also provided, but calculated with the 
same acceptance values as the 7 TeV ATLAS study, and the diffractive to inelastic cross 
section ratio. In Fig. [HI R ss is shown as a function of the diffractive contribution to the 
total inelastic cross section. In the damped scenario, the generator R ss now sits upon the 
central data value. 

In Fig. [13 charged multiplicity distributions at y/s = 7 TeV are shown, giving the 
breakdown of non-, single- and double-diffractive events, (a) over the whole rapidity range 
and (b) for \y\ < 2.5. When the entire phase space is taken into account, the single- and 
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y/s (TeV) 

V V / 


Und (mb) 


crqn (mb) 


<Tr>n (mb) 


Rss (%) 


Crliff /Cinpl (%) 

Hill / Ulcl \ / 




0.90 (pp) 


34.4 


11.7 


6.4 


5.5 


34.5 


Default 


1.96 (pp) 


39.0 


12.5 


7.5 


5.4 


33.8 




7.00 (pp) 


48.5 


13.7 


9.3 


5.1 


32.1 




0.90 (pp) 


36.0 


10.7 


5.8 


5.0 


31.5 


Damped 


1.96 (pp) 


40.9 


11.4 


6.7 


4.8 


30.1 




7.00 (pp) 


50.9 


12.4 


8.1 


4.5 


28.7 



Table 4: Non-, single- and double-diffractive cross sections at different energies for default 
and damped (cr^ x = crg 1 ^ = 65 mb) settings. The percentage of same-side events (using 
the same acceptance values as the 7 TeV study) and ratio of diffractive to inelastic cross 
sections are also given 
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0.06 
» 0.05 
01 0.04 
0.03 
0.02 
0.01 


0.1 0.2 0.3 0.4 0.5 

(a DD + cr SD ) I (a DD + <r SD + cr ND ) 

Figure 14: R ss as a function of diffractive contribution to the total inelastic cross section. 
The dashed curve gives the generator prediction when the ratio is fixed at the default value, 
osd/^dd = 1.53 

double-diffractive events dominate in the lower bins, with the non-diffractive component 
reaching the same level only at N c ^ ~ 30. In (b), just the rapidity cut is enough to 
allow the non-diffractive component to rise above the diffractive one down to much lower 
multiplicities. 

It should be remembered that the non-diffractive cross section enters in the MPI frame- 
work, eq. (126 p . This is most clearly seen in the ratio plots; while the diffractive components 
show a more or less constant drop across the multiplicity range, in (b), the non-diffractive 
component begins to drop in the high- multiplicity tail. This drop will also occur in (a), 
but at higher multiplicities than shown. Any retuning to compensate for the effects of the 
damping must then also take into account the change of slope of the multiplicity distribu- 
tion. 

5.2 Tuning prospects 

We now try a more complete tune to early LHC MB/UE data. Where possible, data has 
been taken from the online HEPDATA database, but certain key datasets are not presently 
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available. For these, instead, data has been read off from the relevant plots available in 
ATLAS publications. In particular, from HEPDATA, the following datasets have been 
taken 

1) ALICE (I77I < 1) charged multiplicity and rapidity distributions at y/s = 
0.90, 2.36 & 7.00 TeV |SQ1 HQ- The charged multiplicity distributions at y/s = 
0.90 & 2.36 TeV are taken from the INEL dataset, but are shown with the zeroth 
bin removed. All others are taken from the INEL>0 dataset, where one track in the 
acceptance region is required to trigger. 

2) ATLAS (\t)\ < 2.5, p± > 500 MeV) charged multiplicity, track p±, mean p± as a 
function of charged multiplicity and rapidity distributions at y/s = 0.90 TeV using 
the INEL>0 dataset [92]. 

The following data has been taken from publications: 

1) As (2) above, but at y/s = 7.00 TeV. Errors for the rapidity distribution are taken, 
but this was not possible for the remaining observables [93J. 

2) ATLAS (\rj\ < 2.5, pj_ > 500 MeV) charged track based underlying event [15]. Charged 
particle number and sum-p^ density distributions are taken with errors for the toward, 
away and transverse regions. A charged track of p± > 1 GeV in the rj acceptance is 
required to trigger an event. 

The results for Tune 2C, introduced in Sec. HI are shown in Fig. (THl compared against 
these datasets. Where available, systematic and statistical errors are shown summed in 
quadrature. At all energies, the multiplicities lie below the data. In the underlying event, 
the charged number density at 7 TeV undershoots in the toward and away regions while at 
both energies there is too little activity in the transverse region. The same plots for Tune 
2M show the same general features. 
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Figure 16: Tune 2C compared against early LHC data (see text). Where available, system- 
atic and statistical errors are shown summed in quadrature as the grey band 



In attempting to improve agreement with these datasets, we stay with the CTEQ6L1 
PDF set, using Tune 2C as a starting point, while turning on the diffractive damping out- 
lined previously. Additionally, we only vary those parameters relating to the MPI framework 
and colour reconnection (specifically the R parameter of Sec. 12.4. 5p . By doing this, we ig- 
nore the possibility of shifting contributions between ISR and MPI; this balance is perhaps 
better determined with observables such asp_i_(Z°) and jet-jet azimuthal correlations. While 
this procedure is unlikely to give us the best tune possible, it should give a good indication 
of whether the generator is capable of describing the broad features of the data. 

The primary tools available, then, are the and E^™ parameters of the MPI frame- 
work, supplemented by £/cm to control the energy scaling (keeping Eq^ fixed at 1800 GeV). 
As an example of the balance that needs to be found between them, Fig. [T7] shows the 
ATLAS and ALICE charged multiplicity distributions and the transverse region charged 
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Figure 17: ATLAS and ALICE multiplicities and the transverse region charged particle 
density when tuned only to the ALICE data (see text) 

particle density when solely tuned to match the ALICE multiplicity data (j/^, = 2.047, 
Ec°p = 1-5 and Eq$[ = 0.19). In order to match the ALICE multiplicities, especially all 
the way out in the tails (made very slightly more difficult by the diffractive damping, as 
shown previously), a lower Ef™ is needed, but this leads to an overshoot in the ATLAS 
multiplicity tails and an underlying event that rises too quickly. 

In Fig. [ISJ results are shown for a set of parameters which give a reasonable description 
to all sets of data, Tune 4C in Tab. [3J In this set, the colour reconnection parameter, R is 
adjusted to match the ATLAS (p±)(N c \ l ) data. This decrease, although preferable in terms 
of describing the data with less non-perturbative corrections, has direct effects on the other 
MB/UE data; it does help to improve the description of the high-multiplicity tails, but the 
extra multiplicity leads to a rise in the transverse charged number density, while there is 
not a similar rise in the J2p± density. The parameter is pushed out such that it in 

fact reduces to the single-Gaussian case. With this setting, the rise in the toward and away 
regions is still slightly too high. While overall, at 7TeV, the multiplicities are slightly on 
the high side, the level of activity in the transverse region in the underlying event is slightly 
too low. As shown in the previous example, getting out to the very high-multiplicity tails 
of the ALICE data is a hard task. 

Fig. shows the charged particle multiplicity at \fs = 1800 GeV and the transverse 
region charged particle density when Tune 4C is run against the Tevatron analyses of the 
previous section. In general Tune 4C gives too much activity, also for the multiplicity 
at y/s = 630 GeV. In both the transverse region (where the J2p± density is in fact well 
described) and overall, the average p±_ per particle sits below the data. This suggests 
a tension between Tevatron and LHC data, or at the very least a new kind of energy 
interpolation between the 900 GeV and 7 TeV data, as suggested by R. Field, based on his 
studies with the Pythia 6 generator [36J. 

6 Summary and outlook 

In this article, modest changes to the parton-shower framework have been introduced, 
including an improved description of dipoles stretching from the final to the initial state 
in the final-state shower and additional azimuthal correlations in the initial-state shower. 
To accompany these changes, studies of the first shower emission in QCD events have been 
made against 2 — > 3 matrix elements, and simple tunes to Tevatron and early LHC data 
have been produced. 
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Figure 18: Tune 4C compared to LHC data 



The azimuthal weighting is primarily based on the expected angular distribution for the 
emission of soft gluons, but it is known that away from this limit it obtains corrections; 
hard emissions are free to "write their own rules". These effects have been studied using 
the a angle, both in comparison to 2 — y 3 matrix elements and Tevatron data. In the 
comparison against data, the width of the pseudorapidity distribution for the softest jet 
is well described, even without the additional weighting, although an improvement in the 
central "dip" is seen. For the a angle, good agreement with 2 — y 3 ME's is noted, but while 
the agreement to data is improved, it is still not as good as hoped for. Further sources of 
asymmetry thus have to be found. 

While an eventual goal would be to have a matching of the first shower emission in QCD 
events to matrix elements, we have begun here with some simple kinematic comparisons. 
Even with such a matching, subsequent emissions would still be handled by the default 
shower, so it is important to understand its behaviour in as much detail as possible. In 
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Figure 19: Charged particle multiplicity at = 1800 GeV and transverse region charged 
particle density at the Tevatron for Tune 4C 

regions where the parton shower is expected to perform well, the description is good both in 
terms of rate and shape, across a wide range of energies. Away from the soft and collinear 
regions, there are non-trivial kinematic effects, but still a reasonable overall agreement, 
generally not worse than could be expected from nontrivial higher orders. 

As shown in Sec. [TJ the inclusion of FSR in the interleaved PS and MPI framework 
has not been without some issues. In particular, the simultaneous description of MB and 
UE data at the Tevatron has been problematic. The changed dipole handling in FSR has 
turned out to be a key factor, and we have performed a simple hand tuning to the Tevatron 
data to show that a combined MB/UE description is possible. This has been done for both 
the CTEQ6L1 and MRST LO** PDF sets. Across all datasets considered, these tunes are 
never significantly worse than the Pro-Q20 and Perugia tunes of Pythia 6.4. These 
tunes are in no way intended to be the final answer; to simplify the task, we have limited 
ourselves to a few key parameters and datasets. It is hoped that what is presented here will 
serve both as a starting point for a more complete global tune to data, and as a general 
guide to some of the principles involved. 

Finally, the first batches of LHC data have been released, including many new MB 
and UE studies. Starting with the CTEQ6L1 Tevatron tune, and limiting ourselves to 
changes in MPI and colour reconnection parameters, supplemented by a slight change in 
diffractive cross sections, a tune that broadly describes the features of much of the data 
has been produced. It is encouraging that this has been possible given the limited scope of 
the tuning. It is hoped that, with a more complete tuning effort, Pythia will be able to 
describe this data more accurately. 

What is shown in this article is a very early look at this new data, and we look forward to 
more data from the LHC experiments. New studies to help separate diffractive contributions 
and their modeling will be welcome. p±(Z°) and jet-jet azimuthal distributions will help 
to separate the contributions from ISR and MPI, while other MB/UE data may help to 
resolve some of the tensions seen so far. It remains an open question as to whether the 
current models will be able to simultaneously describe both Tevatron and LHC data. There 
may well be a region in the overall parameter space that will be able to do this to some 



42 



extent, but we do not rule out differences due to experimental effects or deficiencies in the 
models, e.g. related to the energy dependence of different parameters. 
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